function res = rppca_rfs(rt, gt, pmax, q, indc, j, w)


global xper 

%% INPUT
% Missing values with NaN

% rt          is n by T matrix
% gt          is d by T factor proxies
% pmax        is the upperbound of the number of latent factors
% q           is # of lags used in Newey-West standard errors
% indc        is a 0/1 indicator to include a constant in giglio and xiu
% j           is a number that indicates the observable factor g_t
% w           is the weight to be apply to the means in lettau pelger

%% OUTPUT
% vhatold     estimated latent factors 
% vhat        vhatold imposing positive means and orthogonality
% GammaFM     is Fama and MacBeth estimator
% GammaFM_Crt is Bai and Zhou bias-corrected estimates
% Gammahat    is (d+1) by pmax matrix of risk premia estimates
% avarhatFM   is (d+1) by 1 vector of the avar of risk premia using FM
% avarhatBZ   is (d+1) by 1 vector of the avar of risk premia using BZ
% avarhat     is (d+1) by pmax matrix of the avar of risk premia estimates
% gthat       is d by T by pmax matrix of cleaned factor proxies
% alphahat    is n by pmax matrix of pricing error estimates
% R2G         is d by pmax matrix of time series Rsquare
% R2F         is 1 by pmax matrix of cross-sectional Rsquare

Ir=double(~isnan(rt));
Ig=double(~isnan(gt));

if (sum(Ir(:,1))==0 || Ig(1)==0)
    rt=rt(:, 2:end);
    gt=gt(2:end);
end

if (isnan(gt(1))||isnan(gt(end))) 
   idxnang=find(~isnan(sum(gt,1))); % remove missing values from gt
   gt=gt(idxnang);
   rt=rt(:,idxnang);
end

Ir=double(~isnan(rt));
Ig=double(~isnan(gt));
%% INITIALIZATION

    T  =  size(rt,2);
    n  =  size(rt,1);
    d  =  size(gt,1);
    
    alphahat              =      nan(n,pmax); 
    R2F                   =      nan(1,pmax); %
    R2G                   =      nan(d,pmax); %
    if indc==1
    Gammahat              =      nan(d+1,pmax);
    avarhat               =      nan(d+1,pmax);
    else
    Gammahat              =      nan(d,pmax);
    avarhat               =      nan(d,pmax);
    end     
    
    gthat                 =      nan(d,T,pmax); 
    
    W                     =      nan(pmax,1);
    pvchiW                =      nan(pmax,1);
       
%% ESTIMATION

    rtbar                 =      rt - repmat(mean(rt,2),1,T);% 
    rbar                  =      mean(rt,2);
   
    % Giglio Xiu Method
    gtbar = gt .* Ig;
    for k = 1:d
        observed = (Ig(k,:)==1);
        meanvalue = mean(gtbar(k, observed));
        gtbar(k,observed) = gtbar(k,observed) - meanvalue;
    end        
  %% step 1: PCA
[tempV,tempD]         =      eig((T^(-1))*(rt*rt')+w*(rbar*rbar')); 
[tempeig,tempidx]     =      sort(real(diag(tempD)));
 eigvec                =      tempV(:,tempidx);    
   
                if j==1  
                  MATvarexp = [];
                  for counter=1:pmax
                  varexp=sum(tempeig(end-counter+1:end))/sum(tempeig);
                  MATvarexp = [MATvarexp;  varexp];
                  end
    
                  infoT.fmt = '%10.2f';
                  PClab  = strcat('PC', num2str(repmat(ones(1), pmax, 1)), '-', num2str([1:pmax]'));
                  infoT.cnames = strvcat('Var. Expl.');
                  infoT.rnames = strvcat(' ', PClab);
                end 
 
%% Step 2a: Recover latent factors (vhat) and portfolios' risk exposures (betahat)
 
vhatochk = nan(T, pmax, pmax); 
betaochk = nan(n, pmax, pmax);
vhatchk = nan(T, pmax, pmax); 
betachk = nan(n, pmax, pmax);
equivchck = nan(T, n, pmax);
equivchck2 = nan(T, pmax, pmax);

for phat  = 1 : pmax   % at least one factor up to pmax
     
betahat           =       eigvec(:,(n-phat+1):n);    
vhat              =       betahat'*rt;        
        %need to flip betahat and vhat as the order of the eigenvalues is
        %ascending and not descending (first latent factor is ordered as
        %last, second latend factors is ordered as second to last, etc. notice
        %that this does not impact on Sigmavhat  
        
        betahat               = flip(betahat,2);
        vhat                  = flip(vhat,1);
        
        %make factors orthogonal and rotate loadings accordingly
       
         if phat>1
         [Z,H]=my_orth(vhat');
         vhatold=vhat;
         betahathold=betahat;
         vhat=[];
         vhat=Z;
         rotbeta=inv(H)*betahat';
         betahat=[];
         betahat=rotbeta';
         vhat=vhat';
         
         meanstmp     = mean(vhat, 2);
         signsvtmp    = sign(meanstmp);
         signsvmattmp = signsvtmp * ones(1, size(vhat, 2));
         vhat         = signsvmattmp.*vhat;
         meanstmp     = mean(vhat, 2);
         signsvmattmp = signsvtmp * ones(1, size(betahat, 1));
         betahat      = betahat.*signsvmattmp'; 
         
        end
                
   %% Step 2b: Estimate PC lambdas (gamma) 
        
         % with constant 
         if indc==1
         Gammatilde            =       inv([ones(n,1),betahat]'*[ones(n,1),betahat]) * [ones(n,1),betahat]'*rbar ;
        % without constant
         else           
         Gammatilde            =       inv([betahat]'*[betahat]) * [betahat]'*rbar ; 
         end   
   %% Step 3: Spanning loadings (eta) and mes. error (what)
    
        vhatwithmean          =  vhat;   
        vhat                  =  vhat - nanmean(vhat, 2); 
        Sigmavhat             =  vhat*vhat'/T;  
        etahat               =  gtbar   *  vhat' * inv(vhat*vhat');
        what                  =  gtbar   -  etahat *  vhat;          
   %% step 4: Factor g prices of risk (gamma_g)    

    if indc==1
    Gammahat(:,phat)      =       [Gammatilde(1);etahat*Gammatilde(2:end)]; 
    else
    Gammahat(:,phat)      =       etahat*Gammatilde(1:end);  %added by F.N estimates of the risk premium
    end 
        
 %% step 5: Goodness of Fit
    
        Miota                 =        eye(n) - ones(n,1)*ones(1,n)/n;     
        R2F(1,phat) =  rbar'*Miota*betahat*inv(betahat'*Miota*betahat)*betahat'*Miota*rbar/(rbar'*Miota*rbar);   
        R2G(:,phat)           =        1./diag(gtbar*gtbar').*diag(etahat*(vhat*vhat')*etahat');        
%% step 6: Newey-West Estimation of Avar
        
        Pi11hat               =        zeros(d*phat,d*phat);
        Pi12hat               =        zeros(d*phat,phat);
        Pi22hat               =        zeros(phat,phat);
        
        for t = 1:T
                    
            Pi11hat           =        Pi11hat  +  vec(what(:,t) * vhat(:,t)')*vec(what(:,t) * vhat(:,t)')'/T;
            Pi12hat           =        Pi12hat  +  vec(what(:,t) * vhat(:,t)')*vhat(:,t)'/T;
            Pi22hat           =        Pi22hat  +  vhat(:,t)     * vhat(:,t)'/T;
            
            for s = 1:min(t-1,q) 
                
                Pi11hat       =        Pi11hat  + 1/T*(1-s/(q+1))* (vec(what(:,t) * vhat(:,t)')*vec(what(:,t-s) * vhat(:,t-s)')'+vec(what(:,t-s) * vhat(:,t-s)')*vec(what(:,t) * vhat(:,t)')');
                Pi12hat       =        Pi12hat  + 1/T*(1-s/(q+1))* (vec(what(:,t) * vhat(:,t)')*vhat(:,t-s)' + vec(what(:,t-s) * vhat(:,t-s)')*vhat(:,t)');
                Pi22hat       =        Pi22hat  + 1/T*(1-s/(q+1))* (vhat(:,t)     * vhat(:,t-s)' + vhat(:,t-s) * vhat(:,t)' );            
            end        
        end
        
        % Cross-sectional pricing errors using PCs
        if indc==1
        alphahat(:,phat)   =  rbar  - [ones(n,1),betahat]*Gammatilde; % cross-sectional pricing error using PCs
        else
        alphahat(:,phat)   =  rbar  - betahat*Gammatilde;
        end
        
                     
        if indc==1 
        avarhat(:,phat)   = [(1 - mean(betahat,1)*inv(betahat'*betahat/n)*mean(betahat,1)')^(-1) *var(alphahat(:,phat))/n;...
            diag(kron(Gammatilde(2:end)'*inv(Sigmavhat),eye(d))*Pi11hat*kron(inv(Sigmavhat)*Gammatilde(2:end),eye(d))/T + ...
             kron(Gammatilde(2:end)'*inv(Sigmavhat),eye(d))*Pi12hat*etahat'/T + (kron(Gammatilde(2:end)'*inv(Sigmavhat),eye(d))*Pi12hat*etahat')'/T + ...
             etahat*Pi22hat*etahat'/T) + ...
             diag(var(alphahat(:,phat))*etahat*inv(betahat'*betahat/n - mean(betahat,1)'*mean(betahat,1))*etahat')/n];
       else
       avarhat(:,phat)   = diag(kron(Gammatilde(1:end)'*inv(Sigmavhat),eye(d))*Pi11hat*kron(inv(Sigmavhat)*Gammatilde(1:end),eye(d))/T + ...
       kron(Gammatilde(1:end)'*inv(Sigmavhat),eye(d))*Pi12hat*etahat'/T + (kron(Gammatilde(1:end)'*inv(Sigmavhat),eye(d))*Pi12hat*etahat')'/T + ...
       etahat*Pi22hat*etahat'/T);% + ...
       %diag(var(alphahat(:,phat))*etahat*inv(betahat'*betahat/n - mean(betahat,1)'*mean(betahat,1))*etahat')/n;   
       end
       
%% step 7: P-values gamma_g
       gammat               = Gammahat(:, phat)./sqrt(avarhat(:, phat)); % tstat
       Gammahatp(:, phat)   = f_pval2S(gammat, n-(1+indc))';             % pval
       Gammahatse(:, phat)  = sqrt(avarhat(:, phat));                    % se
       
%% step 8: recovery of gt/cleaned factor (eta x vhat)       
        gthattmp          = etahat*vhat;
        gthat(:, :, phat) = gthattmp; 
 
%% step 9: Wald test for the null that the factor g is weak 
 if d<2
 W(phat) = T*(etahat*(inv(inv(Sigmavhat)*Pi11hat*inv(Sigmavhat)))*etahat');
 pvchiW(phat) = chi2cdf(W(phat),phat,'upper') ;
 else
 W(phat) = NaN;
 pvchiW(phat) = NaN;
 end
    
%% Lettau step (only for j=1 as r and vhat would change only if there are missing obs in gt)
    tempeig_s=sort(tempeig,'descend');
    [MaxSR, weights, iSR, RMS, sigmae, maxf]=lettau_rp(phat, rt', vhatwithmean', tempeig_s); 
    MaxSRmat(:,phat)        = MaxSR;
    RMSmat(:,phat)          = RMS;
    sigmaemat(:,phat)       = sigmae;
    maxfmat(:,phat)         = maxf;
    weightsmat(phat,1:phat) = weights';
    iSRmat(phat,1:phat)     = iSR';
end 
 %% step 11: choose number of latent factors   
    phat        = zeros(3,1);
    ppmax       = 20;
    obj         = tempeig(end:-1:(end-ppmax+1)) + 0.25*(1:ppmax)'* (log(n)+log(T))*(n^(-0.5)+T^(-0.5))*median(tempeig(end:-1:(end-ppmax+1)));
    phat(1)     = find(obj==min(obj))-1; 
    obj2        = tempeig(end:-1:(end-ppmax+1)) + 0.5*(1:ppmax)'* (log(n)+log(T))*(n^(-0.5)+T^(-0.5))*median(tempeig(end:-1:(end-ppmax+1)));
    phat(2)     = find(obj2==min(obj2))-1; 
    obj3        = tempeig(end:-1:(end-ppmax+1)) + 0.75*(1:ppmax)'* (log(n)+log(T))*(n^(-0.5)+T^(-0.5))*median(tempeig(end:-1:(end-ppmax+1)));
    phat(3)     = find(obj3==min(obj3))-1;    
%% OUTPUT

%% latent factors, exposures (both rotated and original), rotating matrix H,sorted eigenvalues,max SR, SDF weights,indivuald SR
res.vhat           = vhatwithmean';
res.betahat        = betahat;   
res.vhatold        = vhatold';
res.betahatold     = betahathold;
res.H              = H;
res.eigenvalues    = sort(tempeig,'descend');
res.MaxSR          = MaxSRmat;
res.RMS            = RMSmat;
res.sigmae         = sigmaemat;
res.weights        = weightsmat;
res.iSR            = iSRmat;
res.maxnf          = maxfmat;
%% cleaned factor, loadings and r2 of spanning regressions
res.etahat         = etahat;      % spanning loadings gt
res.gthat          = gthat;       % gt(c) spanned  factor
res.R2G            = R2G;
%% candidate factor (demeaned) with NaN excludec
res.gtbar          = gtbar;
%% price of risk (gamma_g)
res.Gammahat       = Gammahat;    % gamma_g (coef)
res.Gammahatp      = Gammahatp;   % gamma_g (pval)
res.Gammahatse     = Gammahatse;  % gamma_g (pval)
res.avarhat        = avarhat;     % gamma_g (variance)
%% price of risk latent factors (gammatilde)
res.gammatilde     = Gammatilde;
%% latent factors pricing errors (alphahat), r2 
res.alphahat      = alphahat;  % second pass pricing errors using PCs 
res.R2F           = R2F; 
%% Wald test (statists, pvalue) factor g is weak
res.wstat          = W;
res.wtestp         = pvchiW;
%% Suggested number of factors
res.pmax           = phat;
    
  